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We investigate the fluidization threshold of three-dimensional cohesive granulates under 
hydrodynamic shear forces exerted by a creeping flow. A continuum model of flow through 
porous media provides an analytical expression for the average drag force on a single 
grain. The balance equation for the forces and a force propagation model are then used 
to investigate the effects of porosity and packing structure on the stability of the pile. 
We obtain a closed-form expression for the fluidization threshold of a regular packing of 
mono-disperse frictionless cohesive spherical grains in a planar fracture. The compound 
effect of structural (packing orientation and porosity) and dynamical properties of the 
system on its stability is quantified. 
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1. Introduction 

Granulates are a large collection of macroscopic solid grains. Dry and wet granulates 



2000 



are vital in a large variety of industries, ranging from pharmaceutical to mining (Duran 
pp. 4-10), from construction to agricultural ( Duran ||2000 , pp. 10). They also play 



an important role in many geological processes (" Duran [2000 par. 1.3), such as land and 



mudslides (Dikau fc Brunsden |[l996 ), debris flows (Iverson et al. 2000), erosion, dune 



formation ( Bagnold [1941 ), that shape planets' morphology including, but not limited to. 



Earth ( Hansen et al. |20"TT ). Despite their ubiquity in natural systems, granular materials 
and their behaviour have received considerable attention in the physics community only 
in the last two decades due to their peculiar properties, very different from those com- 
monly associated with either solids, liquids, or gases (Jaeger, Nagel & Behringer ||1996 



|Hornbaker et al. |1997[ |Schiffer |2005[ [Novak, Samadani fc Kudrolli |2005[ and references 
therein) . 

Most theoretical and experimental studies focus on dry granulates and their collec- 



tive behaviour including pattern formation (e.g., Umbanhowar, Melo fc Swinney ||1996 



Hansen et al. 20111), angle of stability /repose (e.g., Jaeger, Liu fc Nagel 1989 Alonso 



fc Herrmann ||1996 ), avalanches dynamics (e .g.,|Jop, Forterre fc Pouliquen ||2QQ6 2007)^ 



and granular flows (e.g., [Rajchenbach ||2003[ |MiDi ||2004^ 

However, as every child knows, adding even a small quantity of liquid to a sandpile 
dramatically changes its properties. The cohesion between grains, due to capillary forces, 
greatly enhances the stability of wet samples. The collective behaviour of wet grains has 



only recently begun to be explored (Jaeger, Nagel fc Behringer 1996 Herminghaus 



2005 and references therein) . The effects of humidity and capillary attractive forces due 



to liquid bridges has a huge impact for practical applications as grains are often not dry 
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in many industrial processes, e.g. when granulates are stored in open air and atmospheric 



humidity plays a role ( Fray fc Marone 200 2), and natural systems, e.g. soil (Iverson 



2000 ). A number of studies have focused on the dynamics of wet granular avalanches (e.g., 
Teg zes, Vicsek & Schiffer '2003) and on the effect of humidity (Fraysse, Thome & Petit 
p999. Fray fc Marone 2002^ Restagno et al ||2002D and capillary forces ([Albert et al 



19971 iHornbaker et al ||l997| [Restagno, Bouquet fc Charlaix ||2004[ [Novak, Samadani fc| 



Kudrolli [2005 ) on the static properties of granulates. A number of different models have 
been proposed to study the geometric stability of wet piles, including Mohr- Coulomb 
continuum ( Mason et al. [[1999 ), liquid-bridge ( Novak, Samadani fc Kudrolli |[20Q5 |) and 
response function ( Ebrahimnazhad-Rahbari et al [[2009 ) models. 

While invaluable in shedding light into the properties of wet granulates, natural sys- 
tems often include a number of additional forcing factors that might significantly affect 
the stability of granulates in the environment. In geological systems, instability is trig- 
gered by a combination of body forces (e.g. gravity) and hydrodynamic shearing due to 
the creeping motion of a fluid through a (often saturated) cohesive granulate. This is 
especially true in processes such as cliffs' instability and landslides, where a fluidization 
event is triggered after heavy rainfall, sediment transport in submerged environments 
(e.g. seafloor transport), fluidization of fines in fractures during pumping operations or 
oil recovery, just to mention a few. Therefore, the understanding of how hydrodynamic 
forces affect the stability of granular matter is of utmost importance to better quantify 
the processes that trigger a fluidization event in submerged granular cohesive systems. 
Moreover, structural properties of cohesive piles, such as porosity, can have a dramatic 
impact on their stability ( Iverson et al. [[2"000[ ). 

In the present work we use a multi-scale framework to quantify the effect of hydrody- 
namic forces on the onset of instability (i.e. fluidization) of a pile constituted of cohesive 
(e.g. oil-wet) grains in a planar fracture. On one hand, a continuum-scale model is used 
to calculate the average hydrodynamic drag exerted on the grains; on the other, a pore- 
scale network model allows us to calculate the force chains propagating through the pile 
and to analytically identify the location of failure and the fuidization threshold of regular 
packing structures for mono-disperse frictionless spherical grains. These are determined 
as a function of structural (e.g. porosity, grain contacts distribution) and dynamical (e.g. 
capillary and hydrodynamic forces) properties of the system. For simplicity, gravity is 
here neglected. Generalisation to include gravity effects is straightforward. 

In section [2] we introduce a continuum-scale Darcy-Brinkmann model for the flow and 
the hydrodynamic force exerted from the creeping fluid on the granulate (Battiato, Ban^ 
daru fc Tartakovsky 2010 Battiato 201 1|. Such model treats the wet granulate as porous 
medium. Therefore, the average drag force exerted on each layer of the pile can be cal- 
culated from continuum-scale equations. Section [3] discusses a pore-scale network model 
for the force propagation through the pile. Location of failure, maximum load and flu- 
idization threshold are analytically derived. Finally, the stability criterion is formulated 
in terms of the capillary number, that represents the relative strength between desta- 
bilising hydrodynamic shear and stabilising capillary (cohesive) forces, and the packing 
orientation relative to the average flow direction. The main results and conclusions are 
summarised in section [4j 



2. Continuum-scale Model of Flow and Drag 

We consider a fully developed incompressible fluid flow between two infinite parallel 
plates separated by the distance of H-\-2L (Fig.[T]). The bottom part of the flow domain, 
—H < ^ < 0, is occupied by a regular packing of (oil-) wet mono-disperse rigid frictionless 




Figure 1. Schematic of the domain on the right, and shape of the average velocity profile 

through the channel, on the left. 



spheres of radius R. The flow is driven by an externally imposed (mean) constant pressure 
gradient dp/dx < 0. Therefore, each spherical grain is subject to a drag force due to 
hydrodynamic stresses and attractive capillary bridge forces. We determine the stability 
regime of this wet three-dimensional sphere packing as a function of hydrodynamic drag 
and packing orientation and structure. 



2.1. Geometry and Packing 

In the present work, we focus on three-dimensional isostatic mono-disperse sphere packing 
obtained by expanding a cubic centred close packing, CCP, so to eliminate interlayer con- 
tacts. Such expanded packing configuration will be referred to as CEP ( Cubic- Expanded- 
Packing). Figures [2| a) and (b) show the side and top view of the structure of the first 
two layers (AB) of spheres. Figure [2|c) shows the top view of a full CEP structure. 

Let us consider a regular pile of two layers of mono-disperse spheres with radius R. 
Let 2R < £ < 2^/3R, be the pitch in the x — z (horizontal) plane, i.e. the distance 
between the centers of spheres belonging to the same layer, and h the_pitch in the y — z 
(vertical) plane, i.e. the distance between two adjacent layers (Fig. |2). When £ = 2R, 
h = 2^/6^/3, the regular close packing is recovered (CCP) and the structure obtained 
by connecting the centers of the four spheres corresponds to a regular tetrahedron (see 
Fig. |2|a)). If £ = 2\^R^ h = and the four spheres lie on the same level. 

Let us define dimensionless quantities 

h , £ 2R ^ , 

h=^, i=j^, e=-. (2.1) 

The dimensionless interlayer and intralayer distances h and £ are related as follows 

h = e\ , €<£< VSe, (2.2) 



and the porosity (j) amounts to 



/. Battiato and J. Vollmer 




Layer A 

Layer B 

Layer C 



Figure 2. Side (a) and top view (b) of the structure of a two-layer 3-dimensional regular pile 
of spheres of diameter e and intralayer distance e < i < y/Se. The distance between the centers 
of two adjacent layers of spheres is h. (c) Top view of a cubic stacking showing the connections 
between spheres' centers belonging to the same layer: from the first two layers of spheres (AB), 
a CEP arrangement can be obtained if every third layer (i.e. ABC ABC.) is the same. 
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€<£< V3e. (2.3) 

When £ = e the close packing is recovered and ^ ^ 0c ^ 0.26, corresponding to a closed 
packing fraction of spheres Sc = 7r/3V^ ^ 0.74. 

Let N -\- 1 he the total number of layers in the pile and n = {!,••• , N} denote the 
layer number. The bottom layer of grains, sitting on the channel's lower boundary, is 
immobile and is called a wall. The layer sitting immediately on the wall has n = 1. The 
layer enumeration continues moving up to the top- most layer in the pile where n = N. 
The number of mobile layers. A/", and the hight of the pile, are related through 
H = Nh + 2R or in terms of dimensionless quantities 

N='-^. (2.4) 

The (continuum-scale) limit e ^ corresponds to N{e) ~ (1/e — 1)/ ^/Y^^]JJe)^J3 oo 
where £/e G (1, V^). 




2.2. Fluid flow model 

Following ( Battiato, Bandaru fc Tartakovsky |2Q1Q Battiato |2Q11 ), we treat the sphere- 
packed region, i.e. the region occupied by spheres, as a porous medium with porosity (j) 
and permeability K. A number of formulas are available to relate permeability to porosity 
for both dilute and non-dilute systems of mono-disperse spheres. In the dilute limit where 
velocity fields of the spheres do not interfere with each other, the permeability is given 
by the well-known Stokes result 



9(1 



(2.5) 



( Torquato 



(1974) and 



2000, pp.505, eq. (19.103)). For nondilute concentration of spheres, Howells 
Hinch I (p^TTl found 



1 



V2 



(1 



a/2 



135 

+ -^(1 - (jy) ln(l - (^) + 16.456(1 - (^) + o(l 



1 -1 



(2.6) 



( Torquato ||2000| p.508, eq. (19.119)). The fiow in this region, y e (-i^,0), can be 
described by the Brinkman equation for the horizontal component of the intrinsic average 
velocity u ( [Auriault ||2009[ ), 



d^2 



K 



dp 
dx 



0, y€{-H,Q), 



(2.7) 



where dp/dx is a mean constant pressure gradient, /i is the fiuid's dynamic viscosity, 
and jie is its "effective" viscosity that accounts for the slip at the spheres walls. Since 
we are concerned with the fiuidization threshold of the sphere packing, initially at rest, 
permeability is constant until the fiuidization threshold is reached. This allows us to 
decouple an analysis of the fiow from that of the granulate dynamics. In the rest of the 
fiow domain, y = [—21/, 0], we use either Navier-Stokes or Reynolds equations to describe 
fully developed fiow respectively in either laminar (7 = 0) or turbulent (7 = 1) regimes 
( P^^[2QQQ1 Eq. 7.8), 

d^^i d{u^v^) dp 
d^^ d^ dx 



p 



0, ^G(0,2L), 



(2.8) 
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where p is the fluid density, and u{y) is the horizontal component of flow velocity u('U, v). 
In the laminar regime, u is the actual velocity and = 0. In the turbulent regime, 
the actual velocity is decomposed into a mean velocity u and velocity fluctuations u' 
and v' about their respective means, {u'v') denotes the Reynolds stress. Fully-developed 
turbulent channel flow has velocity statistics that depend on y only. 

In both flow regimes, the no-slip condition requires zero velocity at y = —H and 
^ = 21/, and the continuity of velocity and shear stress is prescribed at the interface, 
^ = 0, between the free and filtration flows (Vafai fc Kim ||199Qi ): 



u{-H)=0, u{2L)=0, u{0-) =u{0-) = U, pe 



du 
dy 



y=0- 



du 
dy 



(2.9) 



=0+ 



where U is an unknown matching velocity at the interface between channel flow and 
porous medium. 

Let us introduce a characteristic Darcy velocity q = —{H'^/jii)dp/dx and define dimen- 
sionless quantities 



u = 



y 



y_ 



S = 



M = 



Me 



Da = 



K 
IP' 



Rep = 



pHq 

1 



(2.10) 



where Da is the Darcy number (e.g., Niel, Junqueira & Lage 1996, Eq. 6) and Re^ is the 
Reynolds number for porous flow. Then the flow equations (2.7)-([2^ can be rewritten 
in a dimensionless form 

0, 



M 
d^ 



dy^ 
u 

2 ~ 



U 

Da 



1 



7Rep 



d{u'v') 



dy'^ dy 
subject to the boundary conditions 



1 = 0, 



ye (-1,0), 
ye (0,2(5), 



(2.11) 
(2.12) 



u(-l) = 0, w(2(5) = 0, u(O-) = w(0+) = [/, 



du 
dy 



(^ = 0+), (2.13) 



where U = U/q is the unknown dimensionless flow velocity at the interface between the 
free and porous flows. 



The dimensionless version of (2.6) allows one to express the Darcy number Da in terms 
of porosity (j) and the relative diameter of the spheres e. 



Da 



(2.14) 



where 



1 



18(1 - (/.) 



V2 



64 



iln(l 



16.456(1 



+ o{i-4>) 



(2. 



The limit of Da Dac with ^ ^ Dac ~ 1.08e^ and e < 1 (i.e. more than one 
layer of spheres) corresponds to the diminishing flow through the granular pile due to 
decreasing permeability, i.e., due to very dense packing of spheres. In this limit, the 
Brinkman correction term Md'^u/dy'^ in (2.11) becomes negligible and (2.11) reduces 
to Darcy 's law. The limit of Da oo corresponds to free flow and (2.11) reduces to 
the Stokes equation. According to ( |2.14 ) and its graphical representation in Figure |3j 
Da <C 1 for 6 < 1 regardless of the magnitude of (j). This implies that crowding effects 
cannot be neglected for multiple layers of packed spheres. Even in high porosity packing 
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Figure 3. Darcy number Da as a function of porosity and the aspect ratio e (Eq. (2.14)) 



((/) ^ 0.8 — 0.9), common values of e ~ 10 ^ — 10 ^ (e.g. fracture aperture of 100/im and 
clay particle size less than 2/im) give rise to Da ~ 10~^ — 10~^. 

2.3. Hydrodynamic drag on spheres in non-uniform flow 

The drag force per horizontal unit area exerted by the fluid on any cross-section y = const 
is given by the x?/-component of the stress tensor cFxy{y)^ 

du 



(2.16) 



where = /ig for y G [— iY, 0], and = /i for ^ G [0,2L]. The drag force distribution 
along the pile, i.e. the horizontal component of the force D = [^(^),0,0], is obtained as 

-d^i 



(2.17) 



where ^ = 1/i is the number of grains per unit length in the horizontal plane. 

Rewriting ( 2.17| ) in terms of the dimensionless quantities (2.10) and introducing di- 
mensionless shear cr^,,, 



the dimensionless drag ^ is 



H 7/ 

^{y) = m—. 

Ay 



(2.18) 



(2.19) 



The dimensionless drag force F(n) = [F(n),0,0] exerted by the fluid on a sphere be- 



8 /. Battiato and J. Vollmer 

longing to the n-th layer can be obtained by integrating S^{y) over the layer height /i, 



/»-l+(3/i+e)/2 

F(l) = j ^ &{y)dy (2.20a) 

F{i) = / ^{y)dy, z = 2, . • . , iV - 1 (2.206) 

J-l-\-ih 

F{N) = f &{y)dy. (2.20c) 

J-l-\-Nh 

Combining ( |2.20| ) with ( |2.19| ), 

F(l) = mu{3h/2 + e/2 - 1) (2.21a) 

F{i) = m [u{-l + (i + l)h) - u{-l + i/i)] , z = 2, • • • , AT - 1 (2.216) 

F(7V) = [[/ - u{-l + A^/i)] . (2.21c) 

In the following section a pore-scale force network model is used to determine the 
location and the modulus of the maximum load in the pile from the hydrodynamic shear 
exerted by the fluid on each sphere. 



3. Pore-scale Model of Failure 

3.1. Propagation of Forces 

Let F/c be the sum of the external forces acting on grain k (e.g. drag), and gki = Qki^ki 
the force exerted from grain k to grain where hki is a unit vector pointing from grain k 
to grain / and gki is the magnitude of the force; Qki is positive for compressive forces and 
negative otherwise. Whenever there is a stretched capillary bridge between two grains 
k and /, the force exerted by grain k on I is attractive and equal to a constant value 
^ki = fh^ki where > 0. This simplifying assumption that the capillary force is a 
constant, irrespective of grain separation distance, has demonstrated to provide good 



description of collective behaviour of wet granulates ( Ebrahimnazhad-R ahbari et al 
2009| [Fingerle et al ||2008| [Ulrich et ar^2009 ). The force distribution can be uniquely 
determined by solving the following system for the unknowns gki 



yk,l = 1,...,M, 



(3.1) 



with M the total number of grains. There is a unique solution for the force distribu- 



tion (Moukarzel 1998 2004 Ebrahimnazhad-Rahbari et al. 2009) in d dimensions, if 



the packing is isostatic, i.e., if the average number z/ of neighbours per grain equals 2d. 
Note, however, that this solution does not necessarily comply with the constraint that 
contacts break when tensile forces exceed the capillary bridge force. Upon variation of 
parameters a packing yields when this additional requirement is first violated, i.e., the 
packing is stable as long as ^/^^ + > 0, V/c, 

In a regular isostatic 3-dimensional CEP packing of grains bounded by a solid wall 
at the bottom and a free surface at the top each sphere is supported by d = 3 contacts 
from grains further down in the pile. For any couple (/c, /) of grains in contact, the 
unit vector b^^ connecting their centers is aligned to one of the (three) CEP-lattice 
directions. The external force exerted on each grain up in the pile propagates unchanged 
down through the pile to the first mobile layer (n=l) along such directions as discussed 
in Ebrahimnazhad-Rahbari et al. (2009). Therefore, the maximum load is experienced 
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from grains at the bottoni of the pile, i.e. n = 1, and the stabihty threshold is determined 
by the respective bonds carrying the highest load. 

Let F be the maximum force exerted on each grain in the layer n = 1 from the grains 
up in the pile. Let bo,, a = {1,2,3} be the unit vectors along the (lattice) directions 
connecting the center of the supported sphere in the first layer with the centers of the 
supporting spheres of the wall as shown in Fig. [4]^ a). The stabilizing capillary forces 
act along the directions {bi,b2,b3}. If the components of the total force F along such 
directions are less than the capillary forces (assumed constant), the pile is stable. Hence, 
a natural stability criterion can be formulated by decomposing the force F onto the non- 



orthogonal basis {bi,b2,b3} uniquely identified by the structure of the packing (Novak, 
ISamadani fc Kudrolh ||2QQ5D . 



3.1.1. Change of Basis 

Let (o = {ei, e2, 63} be the canonical orthonormal basis of the Euclidean space and 
^ — {bi, b2, b3} a generally non-orthogonal basis with bc^ unit vectors, and a = {1, 2, 3}. 
Let Fq, be the components of the maximum force F in the canonical basis and F' a. its 
components in the basis i.e. 

3 3 

F = ^F,e, = ^F',b,. (3.2) 

Q! = l Q! = l 

The components of F in the two basis are related through a linear transformation A, 

(i^i,i^2,i^3) = A(F\,F'2,n) (3.3) 

with (^1,^2,^3) and {F' \^F' 2^ F' column vectors, and A the matrix of direction cosines 
whose components are defined as A^,^ = cos(b/3, Gq,) = b/3 • Gq,. 

In the ^-coordinate system and for any sphere belonging to the first mobile layer (i.e. 



n = 1), equation (3.1) simplifies to 

F'a = fb a = {1,2, 3}, (3.4) 

and the stability criterion is 

F'o.^fb, a = {1,2, 3}. (3.5) 
Combining (3.3) and (3.5) yields to 

BapFp^h, a, /? = {!, 2, 3}, (3.6) 
where are the components of B := A~^. 

3.1.2. Packing orientation 

The effect of the packing orientation on the pile stability can be explicitly incorporated. 
Without loss of generality, let us consider a counterclockwise rotation of the pile (and 
therefore of the basis ^) around the wall-normal, i.e. ?/-axis (see Figures [4] (a) and (b)). 
Such solid-body rotation is fully described by the rotation matrix ^y{0) defined as 



R,(e) 



COS 6 sin ^ 

1 
— sin 9 cos ( 



(3.7) 



where 6 is the rotation angle. The matrix of direction cosines for the rotated system, A^i, 
is 



(3.8) 
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Figure 4. Top view of the bottom two layers of a regular isostatic packing of mono-disperse 
spheres before (a) and after (b) a counterclockwise rotation of the packing. 



The stability criterion (3.5) now implies 

where Bg^ais are the components of = {^e)~^- 



(3.9) 



3.1.3. Stability criterion: flow in a fracture 

The maximum force F is parallel to the fracture walls. Let the axis e2 of the cartesian 
basis (o be orthogonal to the fracture boundary and ei parallel to F such that F = [F, 0, 0] 
where F = |F|. We take the projection of bi onto the xz-plane parallel to ei as a reference 
configuration for the packing orientation (Fig. [4]^a)). Therefore, the components of the 
basis vectors {bi, b2, bs} in the canonical basis and the matrix of direction cosines A 
are 



e 



-h 






" -Vse/G ' 






-h 




e 


-£/2 






-Vi£/Q - 




-h 


-h 


-h 





-£/2 


£/2 



-V3£/6 
-h 
£/2 



(3.10) 



The matrix of direction cosines A^i after a counterclockwise rotation of angle about 
e2-axis (Fig. [4](b)) has the following components 



(^3^/3) cos (9 -(^3^/6) cos (9 - (^/2) sin i9 

-h -h 
-(V3£/3)sin(9 (V3£/6)sin(9- (£/2)cos(9 



-(V3^/6)cos6> + (^/2) sin^ 
-h 

(V^£/6)sin(9 + (£/2)cos6> 



(3.11) 
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(2V3/3)(e/^)Fcos( 
-(2\/3/3)(e/£)Fcos (6* - 
-(2V3/3)(e/^)Fcos (6*4 



7r/3) 
n/3) 



<fb 
<fb 



(3.12) 



where F is the (as yet unknown) total force exerted by the fluid and the pile on the first 
layer of grains (n = 1). In the following section we determine this force. 

3.2. Shear stress, drag coefficient and maximum load 
As discussed in section 2.3, F(n) is the drag force exerted by the creeping fluid on 
a sphere belonging to layer n. Let F[{n)^ ^2(^) ^si^) be its components along 
the lattice directions bi, b2, bs such that F(n) = F{{n)hi + F2(n)b2 + F3(n)b3. The 
maximum load is experienced at the bottom of the pile where the total force F exerted 
on the grains of the first layer is 



N 



F = ^F(n) = Fi'bi 



F^h2 



(3.13) 



n=l 



where = Y.n=iKi^)^ ^ = {1.2,3}. While equations ( |2.21| ) and ( |3.13[ ) provide a 
general framework for the calculation of the maximum force F in any regular packings 
(e.g. hexagonal), in a CEP configuration F can be readily calculated from (2.19) by 

F = ei/ iMdu = iMUei, (3.14) 



where U is given by (5.3) or (5.5) for laminar or turbulent regime, respectively (see 
Appendix). The dimensional force is therefore F = /HeiU. When Da ^ 0, ~ 5^/D^i/M 
and 



1/2 



which provides a scaling for the maximum force in terms of Da. 

3.3. Fluidization Threshold 



(3.15) 



Combining (3.14) with (3.12), the stability criterion can be written as 



1 - (2V3/3)Ca'^ cos i9 > 
(2V3/3)Ca* cos ((9 - 7r/3) > 
• (2v^/3)Ca* cos((9 + 7r/3) > 



(3.16) 



where 



Ca* 



eMU 



(3.17) 



is a capillary number that incorporates geometrical effects of the porous structure. As- 
suming a toroidal shape of the liquid surface of the capillary bridges, the dimensional 
capillary force = /J^qHfij can be related to the dimensional surface tension 7 = fiqj 
and the contact angle between the wetting liquid (e.g. oil) and the surface of the spheres 
(j) by 



/5 = 27tRj cos (/) 



(3.18) 
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Figure 5. Fluidization threshold of a wet granulate under hydrodynamic shear in terms of 
packing orientation and capillary number Ca"^. 



(Scheel et al. 2008, p. 190) or, in terms of dimensionless quantities, by /5 = 7re7cos( 
Therefore, 



Ca* 



M 



TT COS ( 



-Ca, 



where 



Ca 



7 



(3.19) 



(3.20) 



is the capillary number defined in terms of the interfacial coupling velocity JJ . Solution 
of (3.16) leads the following stability criterion for the CEP as a function of packing 



orientation 9 and capillary number Ca"^, 

Ca* < V3/2, stable V6> 

Ca* > V3, unstable V6> (3.21) 

V^/2 < Ca* < V^, stable for e [arccos(v^/2Ca*), 27r/3 — arccos(v/3/2Ca*)]. 



A graphical representation of (3.21) is provided in Figure [s] The stability of a pile with 
CEP arrangement of its grains is affected by the orientation of the lattice directions 
relative to the average velocity of the flow by a factor of 2 (see Figure [s]): the pile 
orientation determines how the destabilizing hydrodynamic forces decompose along the 
lattice directions and, consequently, how they are balanced by the stabilizing capillary 
forces acting at the contact points. While the height of the pile does not explicitly affect 
its stability, the net shearing velocity U is the parameter that controls the fluidization 
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0, - (^v'Da/M (Eq. (4.3), [Battiato ||2Qll|) that, combined with 
( 3.17| ) and ( |2.14[ ), leads to 

8e^ 



Ca^ = ^VM/((/>). (3.22) 

Therefore, higher porosity piles are less stable than lower porosity ones, as /((/>) is an 
increasing function of (j). This result helps to explain the experimental data of |Iverson et 
ar|(|2000|). 



4. Summary and Conclusions 

By means of a multi-scale framework, we have provided closed-form expressions for the 
fluidization threshold of a regular packing of mono-disperse frictionless cohesive spher- 
ical grains in a planar fracture. The compound effect of structural (e.g. porosity, grain 
contacts distribution, pile orientation) and dynamical (e.g. capillary and hydrodynamic 
forces) properties of the system on its stability is taken into account. The impact of 
packing orientation and porosity is quantified. The orientation of a CEP pile affects its 
stability threshold by a factor of 2 and an increase in porosity renders a pile more un- 
stable. This helps to explain the high impact of porosity on landslides as experimentally 
observed by Iverson et al. ( 2QQQ| . Moreover, while the height of the pile does not explic- 
itly affect its stability, we identify the capillary number, Ca"^ (3.17), as the dimensionless 
parameter that controls the fluidization threshold. Ca"^ is the ratio between destabilizing 
hydrodynamic forces (3.14) and stabilizing cohesive (capillary) forces (3.18). It is de- 
fined in terms of the net shearing velocity t/, the fluid viscosity /i, the surface tension 7, 
and the contact angle (j) between the wetting liquid and the surface of the solid grains. 
For simplicity, gravity was neglected. We stress that, while applied to cohesive capillary 
forces, the method can also be used to model any type of cohesive forces (e.g. Van der 
Waals). The fluidization threshold is then formulated in terms of a dimensionless num- 
ber defined by the same ratio of forces. Generalisation to include gravity effects is also 
straightforward. 

The impact of i) regular packing structures other than cubic (e.g. arrangements derived 
from hexagonal close packing), ii) size and density polidispersity, and iii) randomness in 
grain contacts distribution will be subject of future research. 



5. Appendix 

Solution of the flow problem (2.10)- (2.13) has been derived in ("Battiato, Bandaru & 



Tartakovsky ||2QfQ Battiato"||2Qll ). It is reproduced here for completeness. 



5.1. Flow in the granular pile 

Velocity distribution inside the granular pile is obtained by integrating the Brinkman 
equation (2.11) subject to the first and third boundary conditions ( |2.13 ) 

2x(^) =Da + Cle^^ + C2e-^^ ye [-1,0], (5.1a) 

where A = 1 / VMDa and 

(U - Da)e^ + Da ^ (U - Da)e-^ + Da 



Co 



(5.16) 



The unknown velocity U at the interface y = between the porous and free flows is 
determined by matching the flow in the granular pile (5.1) with the free flow. Velocity 
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profiles for the latter are derived in the following section for laminar and turbulent 
regimes. 



5.2. Flow above the granular pile 



5.2.1. Laminar regime 



boundary conditions (2.13) is 



A solution of the flow equation (2.12) with 7 = subject to the second and third 



u{y) = -'-+[5-^^y + U, yG [0,2,5]. (5.2) 
The interfacial velocity U is obtained from the continuity of the shear stress at = 0, 



the last boundary condition (2.13), as 

^ 1 — sechA 5 tanhA 
U — Da- 



AM P 



/? = 1 



tanhA 
2JAM ' 



(5.3) 



The velocity distribution inside the pile is now uniquely defined by combining (5.1) 



and ( |5.3[ ) . Small Darcy number Da corresponds to a plug- flow regime with uniform veloc- 
ity (i.e. Darcy velocity) almost everywhere in the computational domain; for increasing 
values of Da the velocity profile tends to Poiseuille flow solution (i.e. parabolic profile) 
as shown in Figure 2 from Battiato, Bandaru & Tartakovsky (2010), representing the 
resulting velocity profiles for M 
number Da. 

5.2.2. Turbulent regime 



1, e = 0.001, 5 = 100, and several values of the Darcy 



A number of solutions of (2.12) with 7 = 1 can be found in chapter 7 of Pope (2000). 



Assuming the top of the pile to be hydrodynamically smooth, the dimensionless mean 
velocity u in the viscous sublayer of dimensionless width 5^, obeys the law of the wall 
( P5^[2QQQ1 pp. 270-271), 

u{y)=8y^U, ye [0,8,]. (5.4) 



The continuity of the shear stress at = 0, i.e., the last boundary condition (2.13), yields 
an expression for the dimensionless interfacial velocity, 

8 



U = Da(l - sechA) 



AM 



tanh A. 



(5.5) 



The comparison of (5.3) and (|5.5[) reveals that the interfacial velocities U in the laminar 



and turbulent regimes differ by the factor (3. One can verify that /3 ^ 1 as the Darcy 
number Da 0, i.e. when the permeability and porosity of granular pile become small. 



For such conditions, equation (5.5) represents a good approximation of (5.3). 



REFERENCES 

Albert, R., Albert, L, Hornbaker, D. J., Schiffer, P. & Barabasi, A.-L. 1997 Max- 
imum angle of stability in wet and dry granular media. Phys. Rev E 56(6), R6271-R6274. 

Alonso, J. J. & Herrmann, H. J. 1996 Shape of the tail of a two-dimensional sandpile. 
Phys. Rev. Let. 76(26), 4911-4914. 

AuRlAULT, J.-L. 2009 On the domain of validity of Brinkman's equation. Transp. Porous 
Med. 79, 215-223. 

Bagnold, R. a. 1941 The Physics of Blown Sand and Desert Dunes. Methuen, London, 
UK. 

Battiato, I. 2011 Universal response of carbon nanotube forests to aerodynamic shear. J. 
Fluid Mech. (under review). 



Fluidization of Wet Granulates under Hydrodynamic Shear 



15 



Battiato, I., Bandaru, p. R. & Tartakovsky, D. M. 2010 Elastic response of carbon 
nanotube forests to aerodynamic stresses. Phys. Rev. Lett. 105(14), 144504. 

DiKAU, R. & Brunsden, D. 1996 Landslide Recognition: Identification, Movement and 
Causes. John Wiley & Sons, London, UK. 

DURAN, J. 2000 Sands, Powders, and Grains. Springer- Verlag, New York, USA. 

Ebrahimnazhad-Rahbari, S. H., Vollmer, J., Herminghaus, S. & Brinkmann, M. 2009 
A response function perspective on yielding of wet granular matter. EPL 87, 14002. 

FiNGERLE, A., RoELLER, K., HuANG, K. & Herminghaus, S. 2008 Phase transitions far 
from equilibrium in wet granular matter. New Journal of Physics 10, 053020. 

Fray, K.M. & Marone, C. 2002 Effect of humidity on granular friction at room tempera- 
ture. J. Geophys. Res. 107(B11), 2309. 

Fraysse, N., Thome, H. & Petit, L. 1999 Humidity effects on the stability of a sandpile. 
Eur. Phys. J. B 11, 615-619. 

Hansen, C.J., Bourke, M., Bridges, N.T., Byrne, S., Colon, C, Diniega, S., Dundas, 
C, Herkenhoff, K., McEwen, A., Mellon, M., Portyankina, G. & Thomas, N. 
2011 Seasonal erosion and restoration of mars' northern polar dunes. Science 331, 575- 
578. 

Herminghaus, S. 2005 Dynamics of wet granular matter. Adv. Phys. 54(3), 221-261. 
HiNCH, E.J. 1977 An averaged-equation approach to particle interactions in a fluid suspen- 
sion. J. Fluid Mech. 83, 695-720. 

HORNBAKER, D. J., ALBERT, R., ALBERT, I., BARABASI, A.-L. & SCHIFFER, P. 1997 What 

keeps sandcastles standing? Nature 387, 765. 
HOWELLS, I. D. 1974 Drag due to the motion of a newtonian fluid through a sparse random 

array of small fixed rigid objects. J. Fluid Mech. 64, 449-475. 
IVERSON, R. M. 2000 Landshde triggering by rain infiltration. Water Resour. Res. 36(7), 

1897-1910. 

IVERSON, R. M., Reid, M. E., Iverson, N. R., LaHusen, R. G., Logan, M., Mann, J. E. 
k, Brien, D. L. 2000 Acute sensitivity of landslide rates to initial soil porosity. Science 
290(5491), 513-516. 

Jaeger, H. M., Liu, C. L. & Nagel, R. 1989 Relaxation at the angle of repose. Phys. 
Rev. Let. 62(1), 40-43. 

Jaeger, H. M., Nagel, S. R. & Behringer, R. P. 1996 Granular solids, liquids, and gases. 

Rev. Mod. Phys. 68(4), 1259-1273. 
JOP, P., FORTERRE, Y. & POULIQUEN, O. 2006 A constitutive law for dense granular flows. 

Nature 441, 727. 

JOP, p., FORTERRE, Y. & POULIQUEN, O. 2007 Initiation of granular surface flows in a 
narrow channel. Phys. Fluids 19, 088102. 

Mason, T. G., Levine, A. J., Ertas, D. & Halsey, T. C. 1999 Critical angle of wet 
sandpiles. Phys. Rev. ^60(5), R5044. 

MiDi, GDR 2004 On dense granular flows. Eur. Phys. J. E 14, 341-365. 

Moukarzel, C.F. 1998 Isostatic phase transition and instability in stiff granular materials. 
Phys. Rev. Lett. 81, 1634. 

Moukarzel, C.F. 2004 Carbon nanotubes: science and applications. The Physics of Gran- 
ular Media, Wiley- VCH Verlag. 23-43 pp. 

Niel, D. a., Junqueira, S. L. M. & Lage, J. L. 1996 Forced convection in a fluid-saturated 
porous-medium channel with isothermal or isoflux boundaries. J. Fluid Mech. 322, 201- 
214. 

Novak, S., Samadani, A. & Kudrolli, A. 2005 Maximum angle of stability of a wet gran- 
ular pile. Nature Phys. 1, 50-52. 

Pope, S. B. 2000 Turbulent Flows. Cambridge University Press, New York, USA. 

Rajchenbach, J. 2003 Dense, rapid flowa of inelastic grains under gravity. Phys. Rev. Lett. 
90(14), 144302. 

Restagno, F., Bouquet, L. & Charlaix, E. 2004 Where does a cohesive granular heap 

break? Eur. Phys. J. E 14, 177-183. 
Restagno, F., Ursini, C, Gayvallet, H. & Charlaix, E. 2002 Aging in humid granular 

media. Phys. Rev. E 66, 021304. 



16 /. Battiato and J. Vollmer 

ScHEEL, M., Seemann, R., Brinkmann, M., Michiel, M.di , Sheppard, a., Breidenbach, 
B. & Herminghaus, S. 2008 Morphological clues to wet granular pile stability. Nat 
Mater. 7(189), 189-193. 

Schiffer, p. 2005 A bridge to sandpile stability. Nature Phys. 1, 21-22. 

Tegzes, p., Vicsek, T. & Schiffer, P. 2003 Development of correlations in the dynamics 
of wet granular avalanches. Phys. Rev. E 67, 051303. 

TORQUATO, S. 2000 Random Heterogenous Material - Micro structure and Macroscopic Prop- 
erties. Springer, New York, USA. 

Ulrich, S., Aspelmeier, T., Roeller, K., Fingerle, A., Herminghaus, S. Zippelius, 
A. 2009 Cooling and aggregation in wet granulates. Phys. Rev. Lett. 102, 148002. 

Umbanhowar, P.B., Melo, F. & Swinney, H.L. 1996 Localized excitations in a vertically 
vibrated granular layer. Nature 382, 793-796. 

Vafai, K. & Kim, S. J. 1990 Fluid mechannics of the interface region between a porous 
medium and a fluid layer-an exact solution. Int. J. Heat and Fluid Flow 11(3), 254-256. 



